clear; clc; close all;warning('off')
angstrom = 1e-10;  %[m]
a0 = 5.29177e-11;  %[m]
cmn1 = 0.0299792458e12; %[Hz/cm^-1]
fontsize = 11;
linewidth = 2;
fracY = 0.1;%corresponding to 25 MHz uncertainty of the laser freq relative to 55MHz FWHM
folder = '/Users/MGHu/Documents/KRb REMPI paper_Nat Chem/Paper draft/figure 4/data/';


h1 = figure('position', [0, 500, 700, 600]);
row = 2;
col = 2;
orange = [255 94 8]/norm([255 94 8]);
green = [0, 248, 132]/norm([0, 248, 132]);
blue = [47 108 167]/(47+108+167)*1.9;
lightblue = [0.22, 0.78, 0.99];
red = [102 8 32]/(102+8+32)*1;
lightred = 'r';
gray = [1 1 1]*0.6;

%%%-----plot K2&Rb2 B field theory------
subplot(row,col,1)
fnames = [folder,'B field dependence/', 'MinusFourPlusOneHalf.csv'];
M = csvread(fnames);
Bfield = M(:, 1);
frac_neg4_1over2 = M(:, 2);
alpha_sqr= frac_neg4_1over2;
plot(Bfield, frac_neg4_1over2, '-k', 'LineWidth', 1);
ax = gca;
ax.FontSize = fontsize;
xlabel('B (G)');
ylabel('State admixture probability');
xlim([0, 55]);
ylim([-0.02, 1.02]);
hold on

fnames = [folder,'B field dependence/', 'MinusThreeMinusOneHalf.csv'];
M = csvread(fnames);
Bfield = M(:, 1);
frac_neg3_neg1over2 = M(:, 2);
beta_sqr = frac_neg3_neg1over2;
plot(Bfield, frac_neg3_neg1over2, '--r', 'LineWidth', 1);

fnames = [folder,'B field dependence/', 'MinusTwoMinusThreeHalves.csv'];
M = csvread(fnames);
Bfield = M(:, 1);
frac_neg2_neg3over2 = M(:, 2);
gamma_sqr = frac_neg2_neg3over2;
plot(Bfield, frac_neg2_neg3over2, '-.b', 'LineWidth', 1);
hold off
legend('|\alpha|^2', '|\beta|^2'...
    , '|\gamma|^2', 'Location', 'East');
legend boxoff 

subplot(row,col,2)
PA_B = (alpha_sqr.*beta_sqr + alpha_sqr.*gamma_sqr + beta_sqr.*gamma_sqr);
PS_B = 1-PA_B;
plot(Bfield, PA_B, '-', 'LineWidth', 2, 'Color', gray);
ax = gca;
ax.FontSize = fontsize;
xlabel('B (G)');
ylabel('Probability');
xlim([0, 55]);
ylim([-0.02, 1.02]);
hold on
plot(Bfield, PS_B, '-', 'LineWidth', 2, 'Color', 'k');%[1, 0.5, 0]
hold off
% legend('P_A=|\alpha\beta|^2+|\beta\gamma|^2+|\alpha\gamma|^2', 'P_S=1-P_A', 'Location', 'East');
% legend boxoff 

%%%-----plot K2 B field------
subplot(row,col,3)
fnames = [folder, '20200311_100V_30G_REMPI_scan_5_K2_REMPI.mat'];% K2 J=6

dataxAll = [];
dataK2All = [];
dataK2_bkg_All = [];
dataNcyc = [];
markersize = 8;

load(fnames);
dataxi = xVarList1;                     %[GHz]
dataK2i = xVarCounts_K_2;             %ion counts of K2
dataK2i_bkg = xVarCountsBkgd_K_2;     %bkgnd ion counts of K2
dataNcyci = cycleNumxVar;               %cycle number

dataxAll = dataxi';
dataK2All = dataK2i';
dataK2_bkg_All = dataK2i_bkg';
dataNcyc = dataNcyci';

yplotmin = -0.5;
yplotmax = 30;
Xplot = dataxAll;
Yplot = (dataK2All-dataK2_bkg_All)./dataNcyc;
Yplot_err = sqrt((sqrt(dataK2All+dataK2_bkg_All)./dataNcyc).^2 ...
        +(Yplot*fracY).^2);
if 1
    %---fit data to a Symmetric model----
    xfitdata = Xplot;
    yfitdata = Yplot;
    Aguess = max(yfitdata);
    fun = @(b,x)b(1)*interp1(Bfield, PS_B, x,'spline')+b(2);
    b0 = [Aguess;0];
    weights = @(yhat) 1./((a + b*abs(yhat)).^2);
    [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
    RMSE = sqrt(mean(residual.^2));
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = abs(b - ci(:,1))/2;%%% standard error
    disp('----fit result of K2 vs B----------');
    disp('JK2 = 6');
    disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
    disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
    disp(['RMSE=',num2str(RMSE,3)]);
    curvex1 = Bfield;
    plot(curvex1, fun(b,curvex1), 'Color', 'k',...
        'LineWidth', linewidth,'HandleVisibility','off');
end
    hold on
errorbar(Xplot, Yplot, Yplot_err, 'ok', 'MarkerSize', markersize,...
    'MarkerFaceColor', red,...
    'CapSize', 0, 'LineWidth', 1,'HandleVisibility','off');
plot(Xplot, Yplot, 'ok', 'MarkerSize', markersize,...
    'MarkerFaceColor', red);
xlabel('B (G)');
ylabel('Ion count per cycle');
xlim([0, 55]);
ylim([-0., 15]);

xK2N6 = Xplot';
yK2N6 = Yplot';
yerrK2N6 = Yplot_err';

fnames = [folder, '20200313_100V_30G_REMPI_scan_1_K2_REMPI.mat'];% K2 J=7
dataxAll = [];
dataK2All = [];
dataK2_bkg_All = [];
dataNcyc = [];
load(fnames);
dataxi = xVarList1;                     %[GHz]
dataK2i = xVarCounts_K_2;             %ion counts of K2
dataK2i_bkg = xVarCountsBkgd_K_2;     %bkgnd ion counts of K2
dataNcyci = cycleNumxVar;               %cycle number

dataxAll = dataxi';
dataK2All = dataK2i';
dataK2_bkg_All = dataK2i_bkg';
dataNcyc = dataNcyci';

Xplot = dataxAll;
Yplot = (dataK2All-dataK2_bkg_All)./dataNcyc;
Yplot_err = sqrt((sqrt(dataK2All+dataK2_bkg_All)./dataNcyc).^2 ...
        +(Yplot*fracY).^2);

if 1
    %---fit data to a model----
    xfitdata = Xplot;
    yfitdata = Yplot;
    Aguess = max(yfitdata);
    fun = @(b,x)b(1)*interp1(Bfield, PA_B, x,'spline')+b(2);
    b0 = [Aguess; 0];
    weights = @(yhat) 1./((a + b*abs(yhat)).^2);
    [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
    RMSE = sqrt(mean(residual.^2));
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = abs(b - ci(:,1))/2;%%% standard error
    disp(' ');
    disp('JK2 = 7');
    disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
    disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
    disp(['RMSE=',num2str(RMSE,3)]);
    disp(' ');
    curvex1 = Bfield;
    hold on
    plot(curvex1, fun(b,curvex1), 'Color',gray, 'LineWidth',...
        linewidth,'HandleVisibility','off');
end
errorbar(Xplot, Yplot, Yplot_err, 'sk', 'MarkerSize', markersize,...
   'MarkerFaceColor', blue,...
    'CapSize', 0, 'LineWidth', 1,'HandleVisibility','off');
plot(Xplot, Yplot, 'sk', 'MarkerSize', markersize,...
   'MarkerFaceColor', blue);
ax = gca;
ax.FontSize = fontsize;
hold off
legend('K_2 (N = 6)', 'K_2 (N = 7)', 'Location', 'East');
legend boxoff 

xK2N7 = Xplot';
yK2N7 = Yplot';
yerrK2N7 = Yplot_err';


%%%-----plot Rb2 B field------
subplot(row,col,4)
% fnames = [folder, '20200316_100V_REMPI_scan_Q13_Rb2_REMPI.mat'];% Rb2 J=13
fnames = [folder, '20200317_100V_30G_REMPI_scan_1_Rb2_REMPI.mat'];% Rb2 J=5
dataxAll = [];
dataK2All = [];
dataK2_bkg_All = [];
dataNcyc = [];

load(fnames);
dataxi = xVarList1;                   %[GHz]
dataRb2i = xVarCounts_Rb_2;             %ion counts of Rb2
dataRb2i_bkg = xVarCountsBkgd_Rb_2;     %bkgnd ion counts of Rb2
dataNcyci = cycleNumxVar;               %cycle number

%% ---- the following 6 lines are for averaging J=13 data
fnames = [folder, '20200318_100V_30G_REMPI_scan_1_Rb2_REMPI.mat'];% Rb2 J=13
load(fnames);
dataxi = (dataxi+xVarList1)/2;                     %[GHz]
dataRb2i = dataRb2i+xVarCounts_Rb_2;             %ion counts of Rb2
dataRb2i_bkg = dataRb2i_bkg+xVarCountsBkgd_Rb_2;     %bkgnd ion counts of Rb2
dataNcyci = dataNcyci+cycleNumxVar;               %cycle number

dataxAll = dataxi';
dataRb2All = dataRb2i';
dataRb2_bkg_All = dataRb2i_bkg';
dataNcyc = dataNcyci';

% yplotmin = -0.5;
% yplotmax = 30;
Xplot = dataxAll;
Yplot = (dataRb2All-dataRb2_bkg_All)./dataNcyc;
Yplot_err = sqrt((sqrt(dataRb2All+dataRb2_bkg_All)./dataNcyc).^2 ...
        +(Yplot*fracY).^2);
if 1
    %---fit data to a Asymmetric model----
    xfitdata = Xplot;
    yfitdata = Yplot;
    Aguess = max(yfitdata);
    fun = @(b,x)b(1)*interp1(Bfield, PS_B, x,'spline')+b(2);
    b0 = [Aguess;0];
    weights = @(yhat) 1./((a + b*abs(yhat)).^2);
    [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
    RMSE = sqrt(mean(residual.^2));
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = abs(b - ci(:,1))/2;%%% standard error
    disp('----fit result of Rb2 vs B----------');
    disp('JRb2 = 5');
    disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
    disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
    disp(['RMSE=',num2str(RMSE,3)]);
    curvex1 = Bfield;
%     hold on
    plot(curvex1, fun(b,curvex1), 'Color', 'k', ...
        'LineWidth', linewidth,'HandleVisibility','off');
end
hold on
errorbar(Xplot, Yplot, Yplot_err, 'ok', 'MarkerSize', markersize,...
    'MarkerFaceColor', lightblue,...
    'CapSize', 0, 'LineWidth', 1,'HandleVisibility','off');
plot(Xplot, Yplot, 'ok', 'MarkerSize', markersize,...
    'MarkerFaceColor', lightblue);
xlabel('B (G)');
ylabel('Ion count per cycle');
xlim([0, 55]);
ylim([-0, 24]);
hold off

xRb2N5 = Xplot';
yRb2N5 = Yplot';
yerrRb2N5 = Yplot_err';

% fnames = [folder, '20200316_100V_REMPI_scan_Q14_Rb2_REMPI.mat'];% Rb2 J=14
fnames = [folder, '20200317_100V_30G_REMPI_scan_5_Rb2_REMPI.mat'];% Rb2 J=4
dataxAll = [];
dataK2All = [];
dataK2_bkg_All = [];
dataNcyc = [];
load(fnames);
dataxi = xVarList1;                     %[GHz]
dataK2i = xVarCounts_Rb_2;             %ion counts of K2
dataK2i_bkg = xVarCountsBkgd_Rb_2;     %bkgnd ion counts of K2
dataNcyci = cycleNumxVar;               %cycle number

dataxAll = dataxi';
dataK2All = dataK2i';
dataK2_bkg_All = dataK2i_bkg';
dataNcyc = dataNcyci';

Xplot = dataxAll;
Yplot = (dataK2All-dataK2_bkg_All)./dataNcyc;
Yplot_err = sqrt((sqrt(dataK2All+dataK2_bkg_All)./dataNcyc).^2 ...
        +(Yplot*fracY).^2);
% legend('Rb_2 (J = 13)', 'Rb_2 (J = 14)', 'Location', 'East');
% legend('Rb_2 (N = 5)', 'Rb_2 (N = 4)', 'Location', 'East');
% legend boxoff 
if 1
    %---fit data to a model----
    xfitdata = Xplot;
    yfitdata = Yplot;
    Aguess = max(yfitdata);
    fun = @(b,x)b(1)*interp1(Bfield, PA_B, x,'spline')+b(2);
    b0 = [Aguess;1];
    weights = @(yhat) 1./((a + b*abs(yhat)).^2);
    [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
    RMSE = sqrt(mean(residual.^2));
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = abs(b - ci(:,1))/2;%%% standard error
    disp(' ');
    disp('JRb2 = 4');
    disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
    disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
    disp(['RMSE=',num2str(RMSE,3)]);
    disp(' ');
    curvex1 = Bfield;
    hold on
    plot(curvex1, fun(b,curvex1), 'Color', gray,...
        'LineWidth', linewidth,'HandleVisibility','off');
end
errorbar(Xplot, Yplot, Yplot_err, 'sk', 'MarkerSize', markersize,...
    'MarkerFaceColor', red,...
    'CapSize', 0, 'LineWidth', 1,'HandleVisibility','off');
plot(Xplot, Yplot, 'sk', 'MarkerSize', markersize,...
    'MarkerFaceColor', red);
ax = gca;
ax.FontSize = fontsize;
hold off

xRb2N4 = Xplot';
yRb2N4 = Yplot';
yerrRb2N4 = Yplot_err';



fnames = [folder, '20200316_100V_REMPI_scan_Q13_Rb2_REMPI.mat'];% Rb2 J=13
% fnames = [folder, '20200317_100V_30G_REMPI_scan_1_Rb2_REMPI.mat'];% Rb2 J=5
dataxAll = [];
dataK2All = [];
dataK2_bkg_All = [];
dataNcyc = [];

load(fnames);
dataxi = xVarList1;                   %[GHz]
dataRb2i = xVarCounts_Rb_2;             %ion counts of Rb2
dataRb2i_bkg = xVarCountsBkgd_Rb_2;     %bkgnd ion counts of Rb2
dataNcyci = cycleNumxVar;               %cycle number

%% ---- the following 6 lines are for averaging J=13 data
fnames = [folder, '20200318_100V_30G_REMPI_scan_1_Rb2_REMPI.mat'];% Rb2 J=13
load(fnames);
dataxi = (dataxi+xVarList1)/2;                     %[GHz]
dataRb2i = dataRb2i+xVarCounts_Rb_2;             %ion counts of Rb2
dataRb2i_bkg = dataRb2i_bkg+xVarCountsBkgd_Rb_2;     %bkgnd ion counts of Rb2
dataNcyci = dataNcyci+cycleNumxVar;               %cycle number

dataxAll = dataxi';
dataRb2All = dataRb2i';
dataRb2_bkg_All = dataRb2i_bkg';
dataNcyc = dataNcyci';

% yplotmin = -0.5;
% yplotmax = 30;
Xplot = dataxAll;
Yplot = (dataRb2All-dataRb2_bkg_All)./dataNcyc;
Yplot_err = sqrt((sqrt(dataRb2All+dataRb2_bkg_All)./dataNcyc).^2 ...
        +(Yplot*fracY).^2);
if 1
    %---fit data to a Asymmetric model----
    xfitdata = Xplot;
    yfitdata = Yplot;
    Aguess = max(yfitdata);
    fun = @(b,x)b(1)*interp1(Bfield, PS_B, x,'spline')+b(2);
    b0 = [Aguess;0];
    weights = @(yhat) 1./((a + b*abs(yhat)).^2);
    [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
    RMSE = sqrt(mean(residual.^2));
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = abs(b - ci(:,1))/2;%%% standard error
    disp('----fit result of Rb2 vs B----------');
    disp('JRb2 = 13');
    disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
    disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
    disp(['RMSE=',num2str(RMSE,3)]);
    curvex1 = Bfield;
    hold on
    plot(curvex1, fun(b,curvex1), 'Color', 'k',...
        'LineWidth', linewidth, 'HandleVisibility','off');
end
    hold on
errorbar(Xplot, Yplot, Yplot_err, 'dk', 'MarkerSize', markersize,...
    'MarkerFaceColor', blue, 'CapSize', 0,...
    'LineWidth', 1,'HandleVisibility','off');
plot(Xplot, Yplot, 'dk', 'MarkerSize', markersize,...
    'MarkerFaceColor', blue);
xlabel('B (G)');
ylabel('Ion count per cycle');
xlim([0, 55]);
% ylim([-0.5, 27]);

xRb2N13 = Xplot';
yRb2N13 = Yplot';
yerrRb2N13 = Yplot_err';



fnames = [folder, '20200316_100V_REMPI_scan_Q14_Rb2_REMPI.mat'];% Rb2 J=14
% fnames = [folder, '20200317_100V_30G_REMPI_scan_5_Rb2_REMPI.mat'];% Rb2 J=4
dataxAll = [];
dataK2All = [];
dataK2_bkg_All = [];
dataNcyc = [];
load(fnames);
dataxi = xVarList1;                     %[GHz]
dataK2i = xVarCounts_Rb_2;             %ion counts of K2
dataK2i_bkg = xVarCountsBkgd_Rb_2;     %bkgnd ion counts of K2
dataNcyci = cycleNumxVar;               %cycle number

dataxAll = dataxi';
dataK2All = dataK2i';
dataK2_bkg_All = dataK2i_bkg';
dataNcyc = dataNcyci';

Xplot = dataxAll;
Yplot = (dataK2All-dataK2_bkg_All)./dataNcyc;
Yplot_err = sqrt((sqrt(dataK2All+dataK2_bkg_All)./dataNcyc).^2 ...
        +(Yplot*fracY).^2);

if 1
    %---fit data to a model----
    xfitdata = Xplot;
    yfitdata = Yplot;
    Aguess = max(yfitdata);
    fun = @(b,x)b(1)*interp1(Bfield, PA_B, x,'spline')+b(2);
    b0 = [Aguess;1];
    weights = @(yhat) 1./((a + b*abs(yhat)).^2);
    [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
    RMSE = sqrt(mean(residual.^2));
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = abs(b - ci(:,1))/2;%%% standard error
    disp(' ');
    disp('JRb2 = 14');
    disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
    disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
    disp(['RMSE=',num2str(RMSE,3)]);
    disp(' ');
    curvex1 = Bfield;

    plot(curvex1, fun(b,curvex1), 'Color', gray,...
        'LineWidth', linewidth,'HandleVisibility','off');
end
    hold on
errorbar(Xplot, Yplot, Yplot_err, '^k', 'MarkerSize', markersize,...
    'MarkerFaceColor', lightred,...
    'CapSize', 0, 'LineWidth', 1,'HandleVisibility','off');
plot(Xplot, Yplot, '^k', 'MarkerSize', markersize,...
    'MarkerFaceColor', lightred);
ax = gca;
ax.FontSize = fontsize;
hold off
legend('Rb_2 (N = 5)', 'Rb_2 (N = 4)', 'Rb_2 (N = 13)', 'Rb_2 (N = 14)',...
    'Location', 'East');
legend boxoff 

xRb2N14 = Xplot';
yRb2N14 = Yplot';
yerrRb2N14 = Yplot_err';

% 
% %%%-----plot K2&Rb2 E field theory------
% subplot(row,col,6)
% fnames = [folder,'E field dependence/', 'Efield_MinusFourPlusOneHalfZeroZero.csv'];
% M = csvread(fnames);
% Efield0 = M(:, 1);
% frac_neg4_1over2 = M(:, 2);
% alpha_sqr= frac_neg4_1over2;
% plot(Efield0, frac_neg4_1over2, '-k', 'LineWidth', 2);
% ax = gca;
% ax.FontSize = fontsize;
% xlabel('E (V/cm)');
% ylabel('State admixture probability');
% xlim([0, 400]);
% ylim([-0.02, 1.02]);
% hold on
% 
% fnames = [folder,'E field dependence/', 'Efield_MinusThreeMinusOneHalfZeroZero.csv'];
% M = csvread(fnames);
% Efield0 = M(:, 1);
% frac_neg3_neg1over2 = M(:, 2);
% beta_sqr = frac_neg3_neg1over2;
% plot(Efield0, frac_neg3_neg1over2, '--b', 'LineWidth', 1);
% 
% fnames = [folder,'E field dependence/', 'Efield_MinusTwoMinusThreeHalvesZeroZero.csv'];
% M = csvread(fnames);
% Efield0 = M(:, 1);
% frac_neg2_neg3over2 = M(:, 2);
% gamma_sqr = frac_neg2_neg3over2;
% plot(Efield0, frac_neg2_neg3over2, '-.r', 'LineWidth', 1);
% hold off
% legend('KRb in |-4, 1/2\rangle', 'KRb in |-3, -1/2\rangle'...
%     , 'KRb in |-2, -3/2\rangle', 'Location', 'West');
% legend boxoff 
% PA_E = (alpha_sqr.*beta_sqr + alpha_sqr.*gamma_sqr + beta_sqr.*gamma_sqr);
% PS_E = 1-PA_E;
% 
% %%%-----plot K2 E field------
% subplot(row,col,4)
% Efield = [99-82 990-805 2000-1660]/0.919;   %[V/cm]
% Y_K2_Q7 = [0.65 0.50 0.44];%[ion/cycle]
% Y_K2_Q7_err = sqrt([0.18 0.10 0.09].^2+(Y_K2_Q7*fracY).^2);
% 
% Y_K2_Q6 = [13.21 12.87 10.34];
% Y_K2_Q6_err = sqrt([0.53 0.66 0.54].^2+(Y_K2_Q6*fracY).^2);
% 
% errorbar(Efield, Y_K2_Q6, Y_K2_Q6_err, 'ok', 'MarkerSize', markersize,...
%     'MarkerFaceColor', 'm', 'CapSize', 0, 'LineWidth', 0.7);
% if 1
%     %---fit data to a model----
%     xfitdata = Efield;
%     yfitdata = Y_K2_Q6;
%     Aguess = max(yfitdata);
%     fun = @(b,x)b(1)*interp1(Efield0, PS_E, x,'spline')+0.33+b(2)*0;
%     b0 = [Aguess; 0.33];
%     [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
%     RMSE = sqrt(mean(residual.^2));
%     % Calculate 95% confidence interval of the fitted parameters
%     ci = nlparci(b,residual,'Jacobian',jacobian);   
%     b = b';
%     b_Err = abs(b - ci(:,1))/2;%%% standard error
%     disp('----fit result of K2 vs E----------');
%     disp('JK2=6');
%     disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
%     disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
%     disp(['RMSE=',num2str(RMSE,3)]);
%     curvex1 = Efield0;
%     hold on
%     plot(curvex1, fun(b,curvex1), 'k', 'LineWidth', linewidth,'HandleVisibility','off');
% end
% errorbar(Efield, Y_K2_Q7, Y_K2_Q7_err, 'ok', 'MarkerSize', markersize,...
%     'CapSize', 0, 'LineWidth', 0.7);
% ax = gca;
% ax.FontSize = fontsize;
% xlabel(['E (V/cm)']);
% ylabel('Ion count per cycle');
% xlim([0, 400]);
% ylim([-0.5, 15.5]);
% legend('K_2 (J = 6)', 'K_2 (J = 7)', 'Location', 'West');
% legend boxoff 
% if 1
%     %---fit data to a model----
%     xfitdata = Efield;
%     yfitdata = Y_K2_Q7;
%     Aguess = max(yfitdata);
%     fun = @(b,x)b(1)*interp1(Efield0, PA_E, x,'spline')-0.31+b(2)*0;
%     b0 = [Aguess; -0.31];
%     [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
%     RMSE = sqrt(mean(residual.^2));
%     % Calculate 95% confidence interval of the fitted parameters
%     ci = nlparci(b,residual,'Jacobian',jacobian);   
%     b = b';
%     b_Err = abs(b - ci(:,1))/2;%%% standard error
%     disp(' ');
%     disp('JK2=7');
%     disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
%     disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
%     disp(['RMSE=',num2str(RMSE,3)]);
%     disp(' ');
%     curvex1 = Efield0;
%     hold on
%     plot(curvex1, fun(b,curvex1), 'k', 'LineWidth', linewidth,'HandleVisibility','off');
% end
% 
% %%%-----plot Rb2 E field------
% subplot(row,col,5)
% Efield = [99-82 990-805 2000-1660]/0.919;   %[V/cm]
% 
% Y_Rb2_Q13 = [18.9 20.5 20.8];%Ncyc = 15
% Y_Rb2_Q13_err = sqrt([1.1 1.1 1.2].^2+(Y_Rb2_Q13*fracY).^2);
% 
% Y_Rb2_Q14 = [1.3 1.5 1.1];%Ncyc = 30
% Y_Rb2_Q14_err = sqrt([0.2 0.2 0.2].^2+(Y_Rb2_Q14*fracY).^2);
% 
% errorbar(Efield, Y_Rb2_Q13, Y_Rb2_Q13_err, '^k', 'MarkerSize', markersize,...
%     'MarkerFaceColor', 'c', 'CapSize', 0, 'LineWidth', 0.7);
% if 1
%     %---fit data to a model----
%     xfitdata = Efield;
%     yfitdata = Y_Rb2_Q13;
%     Aguess = max(yfitdata);
%     fun = @(b,x)b(1)*interp1(Efield0, PS_E, x,'spline')+1.57+b(2)*0;
%     b0 = [Aguess; 1.57];
%     [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
%     RMSE = sqrt(mean(residual.^2));
%     % Calculate 95% confidence interval of the fitted parameters
%     ci = nlparci(b,residual,'Jacobian',jacobian);   
%     b = b';
%     b_Err = abs(b - ci(:,1))/2;%%% standard error
%     disp('----fit result of Rb2 vs E----------');
%     disp('JRb2=13');
%     disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
%     disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
%     disp(['RMSE=',num2str(RMSE,3)]);
%     curvex1 = Efield0;
%     hold on
%     plot(curvex1, fun(b,curvex1), 'k', 'LineWidth', linewidth,'HandleVisibility','off');
% end
% 
% errorbar(Efield, Y_Rb2_Q14, Y_Rb2_Q14_err, '^k', 'MarkerSize', markersize,...
%     'CapSize', 0, 'LineWidth', 0.7);
% ax = gca;
% ax.FontSize = fontsize;
% xlabel(['E (V/cm)']);
% ylabel('Ion count per cycle');
% xlim([0, 400]);
% ylim([-0.5, 24.5]);
% legend('Rb_2 (J = 13)', 'Rb_2 (J = 14)', 'Location', 'West');
% legend boxoff 
% if 1
%     %---fit data to a model----
%     xfitdata = Efield;
%     yfitdata = Y_Rb2_Q14;
%     Aguess = max(yfitdata);
%     fun = @(b,x)b(1)*interp1(Efield0, PA_E, x,'spline')+0.67+b(2)*0;
%     b0 = [Aguess;0.67];
%     [b,residual,jacobian,~,MSE,~] = nlinfit(xfitdata,yfitdata,fun,b0);
%     RMSE = sqrt(mean(residual.^2));
%     % Calculate 95% confidence interval of the fitted parameters
%     ci = nlparci(b,residual,'Jacobian',jacobian);   
%     b = b';
%     b_Err = abs(b - ci(:,1))/2;%%% standard error
%     disp(' ');
%     disp('JRb2=14');
%     disp(['A=',num2str(b(1),3), '+-', num2str(b_Err(1),3)]);
%     disp(['y0=',num2str(b(2),3), '+-', num2str(b_Err(2),3)]);
%     disp(['RMSE=',num2str(RMSE,3)]);
%     disp('------------------------------');
%     curvex1 = Efield0;
%     hold on
%     plot(curvex1, fun(b,curvex1), 'k', 'LineWidth', linewidth,'HandleVisibility','off');
% end
set(h1,'Units','Inches');
pos = get(h1,'Position');
set(h1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h1,'plots\BfieldControl','-dpdf','-r0')